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FOREWORD 

This  report  describes  one  phase  of  a program  which  began  in  July  1973  and 
continued  through  June  1975  to  develop  a combined  aero-thermal-stress  computer  code 
capable  of  effeciently  predicting  temperature  and  stress  distributions  in  axisymmetric 
seeker  domes.  The  goal  of  the  program  was  the  determination  of  worst  case 
design/flight  conditions  for  tactical  missile  seeker  domes.  The  program  was  accom- 
plished through  the  Structures  Branch.  Code  4571,  Naval  Weapons  Center,  and 
supported  by  the  Nava)  Air  Systems  Command,  Code  320-B  under  AirTask  A320- 
3 200  008B/5F  3 2-320-206. 
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THIN  is  u special  purpose  program  which  can  be  used  to  obtain  the  stresses 
and  displacements  in  a spherical  shell  acted  upon  by  an  axisynimetric  distribution  ol 
normal  and/or  tangential  pressure,  and  heated  arbitrarily  through  the  shell  thickness  to 
a prescribed  axisynimetric  temperature  distribution.  It  is  assumed  that  the  shell  is 
closed  at  the  apex,  and  supported  along  the  edge  (see  Figure  1)  such  that  the  edge  is 
either: 

1.  Clamped-displacements  (v,  w),  rotation  (X)  all  vanish,  or 


2.  Simply-supported-displacements  (v,  w),  meridional  bending  moment  (M^)  all 


3.  Roller-skate  supported-tangential  displacement  (v),  transverse  shear  (Q^)  and 
meridional  bending  moment  (M^)  all  vanish. 

Further,  it  is  assumed  that  the  material  is  linearly  elastic,  and  characterized  by  material 
parameters  (E.  v,  a)  which  do  not  depend  on  the  temperature.  (This  feature  of  the 
program  could  readily  be  included,  but  was  omitted  here  to  limit  the  length  of  the 
program.) 


FIGURE  1 . Notation  and  Boundary  Conditions. 
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I’HIN  obtains  solutions  of  the  axisymmctric,  thin  shell  equations  appropriate  to 
spherical  shells  (see  Appendix  A).  It  is  based  on  an  algorithm  developed  by  Otter1  and 
teimed  by  him  Dynamic  Relaxation  (D  R.)”.  In  essence,  the  static  problem,  which  is 
governed  by  a system  of  ordinary  differential  equations  and  associated  boundary 
conditions,  is  replaced  by  a dynamic  problem  governed  by  a system  of  partial 
differential  equations  and  the  same  boundary  conditions.  The  success  of  the  method 
depends  on  the  fact  that  the  dynamic  problem  is  relatively  easier  to  solve  numerically 
than  is  the  static  problem. 


The  associated  dynamic  problem  is  governed  by  the  equations  of  motion 
appropriate  to  spherical  shells  in  which  are  included  hypothetical  viscous  damping 
teims.  It  the  damping  parameters  (Kr,  K^)  have  been  chosen  properly,  the  solution  of 
the  equations  ol  motion  should  decay  to  the  required  solution  of  the  equations  of 
equilibrium  in  about  one  to  two  periods  of  the  lowest  mode  of  vibration  of  the  shell 
corresponding  to  the  particular  type  of  support.  (Note  that  rigid  body  motion  must  be 
constrained  in  order  to  conveniently  identify  convergence  of  the  dynamic  solution  to 
the  requited  equilibrium  solution.)  The  use  ol  D. R.  generally  requires  that  the  program 
be  Inst  run  in  an  experimental  mode  with  zero  damping  in  order  to  identify  the 
period  ot  the  lowest  mode  of  vibration.  This  phase  could  be  omitted  if  this 
inlormation  could  conveniently  be  obtained  elsewhere,  i.e,,  analytically.  However,  it  is 
olten  more  easily  obtained  experimentally.  The  required  solution  is  obtained  by 
running  the  program  again  in  the  damped  mode,  and  is  the  limiting  form  of  the  time 
varying  solution  as  time  becomes  large  (about  I 1/2  periods  of  the  lowest  mode  of 
lice  vibration).  In  a qualitative  way,  the  shell  is  first  impulsively  loaded  and  allowed  to 
vibrate  Ireely  (without  damping)  for  a length  of  time  adequate  to  observe  some 
semblance  of  periodic  motion.  This  experiment  is  then  followed  by  a trial  run  in 
which  the  shell  is  again  impulsively  loaded,  but  now  damping  is  included  and  the 
motion  is  observed  for  a length  of  time  adequate  for  the  kinetic  energy  to  be 
dissipated  for  turther  details  ol  D R.,  the  reader  is  referred  to  Rushton.2  In  what 

follows,  the  theory  will  be  implemented  insofar  as  it  applies  to  the  particular  problem 
at  hand. 


IHIN  operates  with  the  thin  shell  equations  made  dimensionless  with  the 
introduction  ol  a suitable  time  scale  and  appropriate  dependent  variables.  These 
equations  (Appendix  A)  are  then  written  in  finite  difference  form  using  the  interlacing 
network  in  both  space  (the  meridional  opening  angle-0)  and  time  that  was  suggested 
by  (lilies-  for  improving  the  accuracy  of  the  finite  difference  representation  of 
derivatives.  Thus,  it  is  observed  from  figures  2 and  3 that  all  the  dependent  variables 
are  not  defined  at  the  same  point  in  either  space  or  time.  In  particular,  the  transverse 


1 Oiler,  .1.  R.  I).  "Computations  for  Pc-s  tressed  Concrete  Pressure  Vessels  Using  Dynamic 
Relaxation,"  Nuclear  Structural  engineering,  / (|%5)  Amsterdam,  pp.  61-75. 

Rushton,  K R.  "Dynamic-Relaxation  Solutions  of  Llastic-Plate  Problems,”  JOURNAL  OF 
STRAIN  ANALYSIS,  Vof,  3,  | (1968),  pp.  23-32. 

_ , , (illics'  D'  ‘Thc  Use  l,f  Interlacing  Nets  for  the  Application  of  Relaxation  Methods  to 
roblems  Involving  Two  Dependent  Variables,"  PROC  ROYAL  SOC.  A (1048),  193,  pp.  407-433. 
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NOMENCLATURE 

SYMBOL  QUANTITIES  EVALUATED 

>«»<•  V *8 
n9;M0.  ms;T 

w,  q 
X,Q 
v.  p 

FIGURE  2.  Spacial  Format. 


QUANTITIES  EVALUATED 
DISPLACEMENTS,  ACCELERATIONS 

VELOCITIES 

DISPLACEMENTS,  ACCELERATIONS 
VELOCITIES 

DISPLACEMENTS,  ACCELERATIONS 
VELOCITIES 


FIGURE  3.  Temporal  Format, 


displacement  (w)  is  defined  on  1 /2-integral  multiple  values  of  the  meridional  coordinate 
interval  (A</>)  and  on  integral  multiple  values  of  the  temporal  coordinate  interval  (At). 
Note  that  the  transverse  velocity  is  defined  on  1 /2-integral  multiple  values  of  the 
temporal  coordinate  interval.  In  comparison,  the  tangential  displacement  (v)  is  defined 
on  integral  multiple  values  of  the  meridional  coordinate  interval,  and  also  on  integral 
multiple  values  of  the  temporal  coordinate  spacing.  In  this  way,  the  representation  of 
derivatives  or  differences  always  involve  central  differences.  On  the  other  hand,  not 
having  quantities  defined  at  points  where  the  boundary  conditions  are  imposed  can  be 
awkward,  However,  this  difficulty  is  overcome  by  the  use  of  dummy  points  which  are 
located  outside  the  regular  interval. 

A typical  solution  is  generated  by  first  assuming  the  shell  to  be  quiescent 
(displacement,  velocity  all  zero)  and  loaded  by  the  given  external  load  and/or  heated 
to  the  prescribed  temperature  distribution.  As  this  condition  can  exist  only  if  there  is 
a non-zero  value  for  the  acceleration,  there  is  motion  away  from  the  initial  configura- 
tion which  can  be  computed  and  continued  in  a step-wise  manner  until  either  the 
period  ot  the  fundamental  mode  of  vibration  is  identified,  or  satisfactory  convergence 
to  the  equilibrium  solution  is  achieved. 

In  the  sections  which  follow,  the  input  requirements  and  the  output  are 
described.  There  then  follows  a description  of  the  program  with  details  on  how  the 
parameters  of  the  program  are  chosen.  Examples  are  presented  of  typical  dynamic 
1 espouse  and  results  given  for  both  constant  transverse  pressure  loading  and  uniform 
temperature  rise.  Finally,  the  accuracy  of  the  program  is  assessed  by  comparing  the 

results  of  1H1N  with  analytical  predictions-again  for  the  case  of  constant  transverse 
pressure, 


INPUT  REQUIREMENTS  AND  THE  OUTPUT  OF  THIN 

In  order  to  use  THIN,  the  following  must  be  provided. 

Material  Data 

v (=PR)  Poisson’s  ratio  (dimensionless) 

a (=ALEA)  Coefficient  of  thermal  expansion  (1/°F) 

(Young’s  Modulus,  E,  is  not  required  as  the  output  stresses  are  presented  in 
dimensionless  form,  i.e„  a H = a00/E). 

Geometric  Data 

0O  (=PO)  Meridional  opening  angle  at  which  boundary  conditions  are 

prescribed  (dimensionless,  radians) 
h (— HB)  Thickness  ratio  (h/R)  (dimensionless) 
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( T lie  thickness,  h,  is  not  required  as  the  output  displacements  are  expressed  in 
dimensionless  form,  i.e.,  w = w/h.) 


A measure  of  the  number  of  segments  (see  Figure  2)  into 
which  the  meridional  length  of  the  shell  is  divided  (N  = 


Control  Data 


Kr  (=C'KR), 
K0  (=CKP) 

NBC 

NV 

MM  AX 


I DEC 


I PR  NT 


Damping  factors  (dimensionless)  (See  Appendix  B) 

An  index  which  determines  (see  Figure  1)  the  boundary 
conditions  (integer) 

An  index  which  determines  the  meridional  station  at  which 
the  tangential  displacement  is  to  be  studied  during  the 
experimental  phase  of  operation  (integer) 

An  index  which  determines  the  number  of  iterations 
through  which  THIN  cycles  in  ei:her  phase  of  operation 
(integer) 

An  index  which  determines  the  interval  between  iterations 
whose  typical  displacements  are  printed  out  during  the 
experimental  phase  of  operation  (integer) 

An  index  which  determines  the  output  required  (integer) 

1PRNT=0  Typical  displacements  every  1DEC  number 
of  iterations  only 

1PRNT=1  Typical  displacements  every  1DEC  number 
of  iterations,  plus  the  solution  after  MMAX 
iterations  expressed  in  dimensionless 
variables 

IPRNT=2  The  solution  after  MMAX  iterations  only 


Mechanical  Loading  and  Temperature  Distribution  Data 

P (=pB)  Dimensionless  tangential  pressure  distribution.  Must  be 

prescribed  at  N-meridional  stations  (see  Figure  2)  even 
though  the  apex  value  is  never  used,  and.  from  symmetry, 
must  be  zero  (subscripted  variable) 

q (=QP8)  Dimensionless  transverse  pressure  distribution.  Must  be 

prescribed  (see  Figure  2)  at  N-meridiona!  stations  even 
though  the  edge  value  is  only  used  for  the  roller-skate 
support  condition  (subscripted  variable) 

TREF  (=TREF)  Reference  temperature  at  which  the  stresses,  strains  are 
assumed  to  vanish  (°R) 
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T (-IMP)  Temperature  distribution.  Temperature  must  be  prescribed 
(see  Figure  4 and  Appendix  C)  at  5 points  through  the 
thickness  of  the  shell,  at  N-meridional  stations  (sub- 
scripted variable,  °R) 

These  data  are  read  into  the  program  from  cards  provided  by  the  user  The 
first  card  contains 


*0  h p Kr  K0 
and  is  read  in  the  open  format.  The  second  card  contains 


and  is  lead  in  the  open  format.  The  third  card  contains 

NBC  N NV  MMAX  1DHC  IPRNT 
nid  is  read  in  the  6110  format.  There  then  follows  N-cards  with  the  number  pr 


0 P 


F = f/h 


pairs: 


.OUTER  SURFACE 


.INNER  SURFACE 


FIGURE  4.  Typical  Temperature  Distribution. 
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for  each  of  N-meridional  stations.  These  cards  are  read  in  the  open  format.  Note  (see 
I igure  2)  that  the  values  of  p,  q are  not  measured  at  the  same  point  in  space  due  to 
the  use  of  the  interlacing  network.  Finally,  there  follows  N-cards  with  the  temperature 
distribution  (see  Figure  4) 

T,  T,  T,  T4  Ts 

through  the  thickness  of  the  shell  for  each  of  N-meridional  stations.  These  cards  are 
also  read  in  the  open  format. 

In  its  present  form,  the  output  of  THIN  always  includes  a record  of  the 
temperature  and  the  applied  pressure  distribution,  and  also  the  temperature  integrals 
1"  . (-)  computed  according  to  Appendix  C.  The  output  can  also  include  (1)  the  time 
response  of  both  the  transverse  displacement  of  the  meridional  station  nearest  the  apex 
(1  = 1)  and  the  tangential  displacement  of  the  meridional  station  1=NV,  and  (2)  the 
solution  after  MMAX  iterations  expressed  in  dimensionless  variables.  The  index  1PRNT 
is  used  to  identify  the  option.  The  displacements  at  1=1  (w),  I=NV(v)  have  been  chosen 
as  they  represent  very  nearly  the  maximum  displacements  expected  in  the  lowest  mode 
of  vibration.  The  index  NV  is  taken  somewhat  greater  than  N/2  as  it  has  been 
observed  that  the  tangential  displacement  reaches  a peak  in  this  neighborhood.  The 
time  response  is  printed  out  so  that  the  user  can  either  estimate  the  period  of  the 
lowest  mode  of  vibration  or  determine  whether  satisfactory  convergence  to  the  static 
solution  has  been  achieved.  The  solution  includes  the  meridional  distribution  of  the 
dimensionless  displacements  and  stress  resultants,  and  the  complete  (meridional  plus 
transverse)  distribution  of  the  dimensionless  meridional  normal  stress  (a^  = o^/E). 
This  stress  component  is  computed  according  to  the  method  outlined  in  Appendix  D, 
and  is  displayed  as  it  is  generally  the  critical  stress  from  a maximum  normal  stress 
point  of  view. 

ANALYSIS  OF  THIN 

The  form  of  the  thin  shell  equations  used  in  THIN  is  that  given  by  Parkus  and 
I'lugge.4  This  form  was  chosen  as  it  was  considered  more  generally  accessible  than 
other  forms.  The  equations  and  their  dimensionless  forms  are  presented  in  Appendix 
A.  THIN  actually  operates  with  the  Finite  difference  form  of  these  equations  using  the 
interlacing  network  illustrated  in  Figure  2. 

The  D.R.  algorithm  can  be  illustrated  by  studying  the  equations  which  govern 
the  transverse  displacement  (w).  _At  the  beginning  of  a typical  iteration,  the  displace- 

3w 

merits  w(I)  and  the  velocities  ^r(I)  are  known  at  all  points  in  space  (1<I<N).  It 
should  be  noted  (see  Figure  3)  that  the  velocities  are  known  at  a time  (corresponding 


4 Fluggc,  W.,  Handbook  of  Engineering  Mechanics,  McGraw-Hill,  1962. 
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to  the  index  M)  which  is  less  than  the  time  (also  corresponding  to  the  index  M)  at 
which  the  displacements  are  known.  For  the  first  iteration  (M=l),  these  displacements 
and  velocities  are  all  taken  equal  to  zero.  Before  proceeding  to  the  next  iteration,  the 
strain  components  e^,  eg  are  first  computed.  For  example,  from  Appendix  A,  we  write 

e0(i)  = w(I)  + [v(I+l ) - v(l)]/A0 

In  an  analogous  manner,  the  rotation  X is  also  computed.  With  the  strain,  rotation 
components  now  known  at  all  points  in  space,  the  in-plane  force  and  moment  stress 
resultants  follow  directly  from  the  constitutive  equations-the  temperature  integrals 
having  been  computed  (see  Appendix  C)  during  the  data  loading  phase.  The  transverse 
stress  resultant  (Q)  is  then  obtained  from  the  equation  of  moment  equilibrium 
expressed  in  the  form 


Q(D  = IM0(I)  - M0(I-  1)]/A0  + [M0(I)  - M0(l)  + M0(I-  1)  - M0(I-1)] 

Note  that  the  use  of  this  equation  is,  strictly  speaking,  not  justified  in  the  context  of 
a dynamic  analysis.  However,  since  the  dynamical  aspect  used  here  is  a computational 
artifice,  the  neglect  of  rotary  inertia  is  unimportant  and  has  no  effect  on  the  final 
results.  With  all  stress  resultants  now  known,  the  solution  corresponding  to  the  given 
displacements  is  complete  and  the  stress  distribution  could  be  calculated  although  it  is 
generally  deferred  until  the  last  (M=MMAX)  iteration.  It  should  be  noted  that  this 
solution  does  not  in  general  satisfy  the  equation  of  translational  equilibrium.  It  is 
precisely  this  fact  that  enables  the  calculation  to  proceed. 

rite  displacement  distribution  for  the  next  iteration  (corresponding  to  the  index 
M+l)  is  obtained  from  the  acceleration  which  is  computed  from  the  translational 
equation  of  motion.  This  requires  the  definition  of  the  acceleration  expressed  in  terms 
of  the  velocities  given  by 

?)“w  3w  3w 

^T(M)  = [-~(M+D-  ^(M)l/Ar 

and  the  definition  of  the  velocity  in  terms  of  the  displacement  given  by 
|^(M+1)=  [w(M+l)~  w(M)l/Ar 

Thus,  with  the  transverse  translational  equation  of  motion  expressed  in  the  form 

^(M+l)(l  + Kr/2)  = |y(M)(  1 - Kr/2)  + Arq(I)-  (N0(l)  + Nfl(l) + 

+ ((0(1+1)  - Q(1))/A0  + Sy-^(Q(1+ 1 ) + Q(I)))) 

1 2R"  ^ 
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the  velocities  ^(I.M+1)  and  the  displacements  w(l,M+l)  follow  directly.  This 

completes  one  cycle  of  the  D.R.  algorithm  save  for  the  establishment  of  dependent 
variables  at  dummy  points  that  are  necessary  for  the  satisfaction  of  the  boundary 
conditions. 


l-ollowing  the  format  given  on  Figure  2,  it  is  apparent  that  boundary 
conditions  can  be  directly  enforced  on  w,  as  they  are  defined  on  0 = 0O_ 

Alternatively,  the  requirement  that  the  tangential  displacement  vanish  at  <p  = <p  leads 
to  the  establishment  of  v at  the  dummy  point  1=N+1  with  the  valise 

vfN+1 ) = -v(N) 

Thus  it  is  actually  the  average  of  two  adjacent  points  which  actually  is  required  to 
vanish.  A similar  expression  is  used  to  enforce  a zero  condition  on  variables  which  are 
not  defined  on  0 = 0(), 


The  clamped  boundary  condition  is  enforced  by  taking 
X(N+1 ) = - X(N) 

It  then  follows  that 

k0(N)  = 4 (X(N+l)  + X(N))cot  0Q  =0 
*0<N)  = (X(N+1)  - X(N))/A0  = -2X(N)/A0 

which  is  sufficient  to  determine  the  edge  values  of  the  moment  stress  resultants. 

The  simply-supported  boundary  condition  is  enforced  by  taking 
M0(N)  = O 

This  requirement  is  used  to  establish  X(N+!)  from  the  constitutive  equation  expressed 
in  the  form 


M0(N)  = (X(N+1)-  X(N))/A0  + -j  (X(N+1)  + X(N))  cot  0„  + ( l+v)  «0R/h  = 0 
Thus,  it  follows  that 


X(N+ 1 ) = (X(N)f  1 - ^ cot  0Q ) - ©A 0)/(  1 + cot  0O ) 


11 
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The  roller-skate  boundary  condition  is  the  only  option  in  THIN  in  which 
wtNTAO,  As  the  roller-skate  is  tree  to  move  in  the  transverse  direction,  the  transverse 
shear  mush  vanish  leading  to  a dummy  value  for  Q given  by 

QtN+1)  =-Q(N) 

With  Q(N+1)  now  known,  the  transverse  displacement  at  the  edge  can  be  computed 
using  the  transverse  translational  equation  of  motion. 


Till  US!-!  OF  THIN  - EXAMPLES 

The  primary  unknowns  in  the  D.R.  agorithm  are  the  damping  factors  Kf,  K0. 
In  order  to  determine  these  factors  (see  Appendix  B),  the  program  must  generally  first 
be  run  in  the  undamped  mode  to  estimate  the  lowest  natural  frequency.  However, 
before  this  can  be  done,  the  user  must  choose  an  appropriate  value  for  N and  estimate 
the  total  number  of  iterations  required  (MMAX)  to  identify  the  period  of  the  lowest 
mode  of  vibration. 

The  index  N is  an  integer  and  is  a measure  of  size  of  the  spacial  interval  A 0.  It 
is  defined  by 


The  choice  of  N is  an  important  one,  since  it  not  only  affects  the  accuracy  of  the 
solution,  but  essentially  establishes  the  running  time  of  the  program  and  hence  the 
cost.  If  N is  large,  the  finite  difference  representation  of  the  spacial  derivatives  is 
relatively  accurate.  However,  the  allowable  time  step  (Ar)  is  limited  by  the  require- 
ments that  the  algorithm  be  numerically  stable  and  the  program  must  be  cycled  many 
times  in  order  that  the  overall  time  interval  be  approximately  1 1/2  periods  of  the 
lowest  mode.  Alternatively,  if  N is  small,  the  finite  difference  representation  becomes 
increasingly  inaccurate,  but  the  running  time  is  substantially  reduced.  In  the  final 
analysis,  the  choice  of  N can  only  be  made  on  the  basis  of  several  trials  with  the 
criteria  being  that  N is  chosen  to  be  the  smallest  number  such  that  there  is  a 
negligible  change  in  the  solution  corresponding  to  a finite  increase  in  N. 

As  a general  rule,  the  spacing  along  the  meridian  RA0  should  be  less  than  the 
characteristic  length  of  the  shell  (^hR).  If  we  chose  the  spacing  A0  such  that 

RA0  =<v/hR/3 


30o 

then,  N =«  and  we  should  obtain  about  ten  data 
v h 

corresponding  to  three  characteristic  lengths. 


points  in  a meridional  distance 


12 


NWC  TP  5785 


The  choice  of  MMAX  depends  on  the  allowable  time  increment  A r.  As  the 
equations  of  motion  governing  the  spherical  shell  are  probably  of  the  parabolic  type, 
we  adopt  the  requirement  established  by  Crandalls  for  beams  that  limits  the  stable 
time  step  to 


A r < W 

f urther,  if  we  assume  that  the  period  of  the  lowest  mode  of  vibration  is  approx- 
imately (r()  =)27r,  and  adopt  the  equality  in  the  above  expression  for  the  stable  time 
step,  a running  time  of  a period  and  a half  corresponds  to 

MMAX  = 6ttN2/05 

for  example,  for  a hemispherical  shell  (0Q  = 90  deg)  and  N=10,  a running  time  of  a 
period  and  a half  corresponds  to  MMAX=764.  With  N,  MMAX  estimated  approx- 
imately. the  program  may  be  run  and  improved  values  obtained  from  an  examination 
of  the  solution  and  the  time  behavior. 

Hxamples  of  undamped  response  are  presented  as  Figures  5 through  8 for 
constant  normal  pressure  loading.  (The  response  to  a constant  temperature  rise  for  the 
same  boundary  conditions  and  shell  geometry  is  not  shown  since  the  curves  are  similar 
to  those  for  constant  pressure  loading.)  Note  that  as  the  curves  have  been  plotted  with 
the  iteration  index  (M)  as  the  abscissa  (and  with  an  iteration  step  size  1DEC=40)  the 
corresponding  dimensionless  time  is  different  for  each  shell  geometry  and  must  be 
computed  according  to  the  relation 

t = MAt  = (A0)“ 

Some  similarity  can  be  seen  by  comparing  the  curves  on  Figures  5 and  6.  It  appears 
that  the  effect  of  the  thickness  parameter  (h)  is  greater  for  the  tangential  displacement 
response,  as  the  transverse  displacement  response  have  roughly  the  same  shape.  One 
should  hasten  to  add,  however,  that  there  is  not  complete  geometric  similarity  between 
the  shells  as  the  tangential  displacement  was  observed  at  0 « 52  deg  for  the  case  N=16 
and  at  0 ^ 42  deg  for  the  ease  N=20.  Nevertheless,  it  appears  that  the  use  of  r as  a 
time  scale  is  justified,  and  that  some  semblance  of  periodic  motion  may  be  observed 
corresponding  to  a dimensionless  period  of  about  6<r<8.  Thus,  a running  time 
(MMAX)  corresponding  to  9 seems  to  be  adequate  to  span  about  1 1/2  periods  of 
the  lowest  mode  of  vibration  for  clamped  boundary  conditions  at  0O  = 90  deg. 


5 Crandall,  S.  H.  Engineering  Analysis,  A Survey  of  Numerical  Procedures,  McGraw-Hill,  1956, 
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DIMENSIONLESS  TRANSVERSE  DISPLACEMENT 
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IGURE  7.  Transverse  Displacement  Versus  Time. 
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FIGURE  8.  Tangential  Displacement  Versus  Time. 

The  cl  foots  of  boundary  condition  and  opening  angle  can  be  seen  by  comparing 
the  responses  shown  on  Figures  5 and  7.  It  is  apparent  that  the  stiffening  effect  of 
lie  clamped  support  at  0O  - 35  deg  substantially  lowers  the  period  of  the  lowest 
mode  ol  vibration,  further,  there  are  only  two  local  maximas  in  the  35-deg  shell 
response  prior  to  crossing  the  axis  in  comparison  to  three  local  maxima  for  the  90-deg 

shell.  1 inis,  a running  time  corresponding  to  r = 8 should  be  adequate  to  span  1 W 
periods  for  a shell  clamped  at  0()  = 35  deg. 

It  can  be  shown  that  the  effect  of  replacing  the  clamped  support  by  a simple 
support  at  the  same  value  of  opening  angle  and  thickness  ratio  is  insignificant  and  does 
not  warrant  any  special  attention.  However,  as  can  be  seen  from  Figure  7.  a dramatic 
el  feet  is  achieved  m the  transverse  displacement  response  if  one  replaces  the  clamped 
support  at  0n  - 90  deg  with  a roller-skate  support.  As  the  tangential  displacement  is 
identically  zero,  the  breathing  mode  only  has  been  excited,  and  this  corresponds  to  a 
dimensionless  period  of  approximately  400  (±20)  X Ar  = 3.73±0.19.  This  result 
compares  well  with  the  prediction  given  by  Kraus6  that 


Kraus,  H.  Thin  Elastic  Shells , John  Wiley  Book  Co.,  1967. 
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different,  based  on 


observation  that  the  frequency  of  the  coupled  mode  was  lower  than  that  of  the 
breathing  mode  Thus,  Kf  K0  are  chosen  equal  and  inversely  proportional  to  the  lowest 
period  of  vibration  in  all  calculations. 

At  this  stage,  one  can  choose  the  output  control  parameter  IPRNT=1  and 

obtain  the  complete  solution  at  M=MMAX,  and  the  damped  displacement  response  to 

use  as  a check  for  convergence.  However,  if  one  has  sufficient  confidence  in  the  choice 
of  damping  coefficient,  1PRNT  can  be  set  equal  to  2 and  only  the  final  solution  will 
be  obtained. 

Examples  of  displacement  and  stress  distributions  are  presented  in  Figures  10 
through  12  for  a shell  of  thickness  ratio  R = 1/8,  clamped  at  (pQ  = 90  deg  and  heated 
uniformly  (Tm  = 0.060).  As  expected,  the  results  indicate  an  edge-effect;  the 

significant  stress  distribution  is  confined  to  a region  near  the  restrained  edge. 

Quantitatively,  the  stresses  seem  to  die  out  in  a distance  of  about  three  characteristic 
lengths  (s/hR)  from  the  clamped  edge.  This  confirms  the  recommendation  made  earlier 
in  choosing  a value  of  N. 


I 
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Accuracy  of  THIN 

I 

The  accuracy  of  THIN  can  be  assessed  by  comparing  its  predictions  with  exact 
solutions  of  the  governing  equations.  However,  since  there  are  very  few  exact  solutions 
available  (which  are  of  any  real  interest),  there  are  a number  of  solutions  whose 
accuracy  can  be  estimated  and  therefore  can  qualify  as  a standard  for  comparison. 

ii 

One  exact  solution  which  is  meaningful  is  the  membrane  state  corresponding  to 
uniform  normal  pressure  and  supported  at  the  edge  so  as  to  restrain  the  tangential 
displacement  and  not  induce  any  bending.  This  solution  can  be  expressed  in  the  form 


l-n_  — — l-j/2_  _ 

w = -y-q  N0,N0=-~2-q  v,Q,M^,M0=O 


This  state  can  be  achieved  with  THIN  by  using  the  roller-skate  support,  and  leads  to 
w = 0. 1440X  10" 2 N0,N0  =-0.1800  X 10' 2 M0,M0  =0.000 

for  all  stations.  The  program  parameters  used  in  this  solution  were 

N = 12  $0  = 90  deg  v=  1/4  fi=l/8  q = -0.384  X 10~2 

Kr,  K0  = 0.0220  MMAX  = 1000 
As  can  be  seen,  the  agreement  is  perfect. 
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CONSTANT  TEMPERATURE  (Tr 
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FIGURE  10.  Transverse  Displacement  Distribution 
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FIGURE  11.  Tangential  Displacement  Distribution, 
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Besides  this  simple  example,  there  are  (to  this  writer’s  knowledge)  no  other 
convenient  analytical  solutions  available.  There  are,  however,  some  series  solutions 
whose  accuracy  should  improve  with  the  number  of  terms  considered.  These  solutions 
tall  into  two  categories:  ( 1 ) asymptotic  solutions  developed  in  powers  of  the  thickness 
parameter  (h)  and  (2)  series  solutions  developed  in  powers  of  the  meridional 
coordinate  (0). 

In  the  first  category,  it  is  a straightforward  task  to  develop  the  solution  for  a 
clamped  hemispherical  shell  acted  upon  by  uniform  normal  pressure.  The  details  of  the 
procedure  are  presented  in  Appendix  E and  only  the  numerical  results  will  be 
discussed  here.  Although  the  convergence  of  the  solution  to  the  exact  solution  has  not 
(to  this  writer’s  knowledge)  been  established,  one  can  at  least  assert  that  the  error  in 
the  solution  should  decrease  as  the  thickness  parameter  (h)  is  decreased.  For  very  small 
h.  we  would  expect  that  the  difference  between  the  (unknown)  exact  solution  and  the 
two-term  asymptotic  solution  presented  here  should  also  be  small.  Thus,  if  the  results 
predicted  by  THIN  compare  favorably  with  those  of  the  asymptotic  solution,  one 
could  conclude  that  the  error  in  the  numerical  process  was  also  small. 

Numerical  results  are  presented  in  Figures  13  through  15  for  the  displacements 
and  circumferential  normal  stress  resultant  obtained  from  THIN  and  from  a two-term, 
unitormly  valid  asymptotic  solution  (UVS)  corresponding  to  thickness  parameters  (fi) 
ot  1/8,  1/16,  0.033.  The  meridional  normal  stress  resultant  (N^/Rpr)  was  not  plotted 
since  both  the  predictions  of  THIN  and  the  UVS  varied  only  insignificantly  from  each 
other  and  from  the  value  0.5000  over  the  entire  range  (O<0<9O  deg).  As  can  be 
seen,  the  agreement  between  the  two  methods  is  generally  satisfactory.  The  disparity 
between  the  two  methods  for  the  transverse  displacement  seems  to  increase  with 
decreasing  values  of  R.  Alternatively,  the  disparity  in  the  tangential  displacement  seems 
not  to  vary  with  R,  and  the  UVS  for  the  circumferential  normal  stress  resultant  is  not 
shown  as  it  differed  imperceptibly  from  the  results  of  THIN.  The  behavior  of  the 
displacement  results  can  be  understood  when  considered  along  with  the  results 
presented  in  Figure  9.  As  can  be  seen  in  Figure  9,  the  amplitude  of  the  oscillation  still 
existing  when  the  iterations  were  terminated  was  of  order  3%.  Further,  the  final  results 
stem  to  depend  on  the  choice  of  the  damping  parameters.  Thus,  although  the  accuracy 
of  the  UVS  should  be  increasing  with  decreasing  value  of  h,  the  accuracy  of  THIN 
would  always  be  expected  to  be  limited  to  the  range  2-3%.  It  is  noteworthy  that  the 
agreement  between  the  results  for  the  stress  resultant  N^  does  not  vary  with  h. 

For  a solution  in  the  power  series  category,  one  can  refer  to  the  example  given 
by  Timoshenko  and  Woinowsky-Krieger7  for  a spherical  cap  loaded  by  uniform  normal 
pressure  and  clamped  at  0Q  = 35  deg.  The  analytical  details  are  described  and  curves 
presented  for  the  distribution  of  meridional  bending  moment  and  circumferential 
normal  stress  resultant  (due  to  bending).  It  is  noted  that  the  numerical  residts  were 
obtained  by  summing  ten  terms  of  the  series  solution.  The  implication  seems  to  be 


7 Timoshenko,  S.  P.  and  S.  Woinowsky-Krieger.  Theory  of  Plates  and  Shells,  2nd  F.di'.ion, 
McGraw-Hill,  1959. 
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TWO-TERM  ASYMPTOTIC  EXPANSION 
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that  the  authors  considered  this  to  be  adequate,  and,  therefore,  these  results  should  be 
useful  as  a standard  of  comparison.  Similar  results  have  been  obtained  using  THIN  and 
are  presented  in  Figure  16  for  the  following  program  parameters: 

N = 8 <pQ  = 35  deg  i>=1/6  h = 1/30  MMAX=  2400 
q = ~ 0. 3 X 1 0"  4 Kr,K0  = 1.77  X At 


No  attempt  has  been  made  to  reproduce  the  curves  from  Timoshenko’s  work  since  the 
figures  there  are  small  and  could  only  be  reproduced  with  difficulty  and  with  dubious 
accuracy.  It  is  therefore  left  to  the  reader  to  confirm  the  observation  that  the 
agreement  is  satisfactory.  It  should  be  noted,  however,  that  the  results  of  THIN  for 
the  meridional  bending  moment  seem  to  underestimate  the  edge  value  by  about  10%. 

The  general  conclusion  from  the  three  examples  presented  above  is  that  THIN 
offers  an  accurate  means  of  obtaining  the  stresses  and  displacements  in  spherical  shells 
with  a variety  of  boundary  conditions  enforced  on  edges  ranging  from  shallow  to 
hemispherical.  Admittedly,  the  loading  conditions  studied  here  (uniform  pressure, 
uniform  temperature  rise)  were  simple.  However,  the  restraint  applied  at  the  edge  in 
each  case  (see  Figure  15,  for  example)  produced  a boundary  layer  (rapidly  varying 
stress  state)  that  was  accurately  predicted  by  THIN.  As  this  boundary  layer  can  be 
regarded  as  an  example  of  a singular  solution  out  of  which  more  complex  stress  states 
can  be  built  up  by  superposition,  one  can  conclude  that  THIN  has  the  capacity  of 
predicting  the  stress  states  due  to  more  complex  loading  conditions  with  an  accuracy 
at  least  comparable  to  that  obtained  above. 
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Appendix  A 

THE  GOVERNING  EQUATIONS 

The  governing  equations  are  taken  from  sections  of  the  Handbook  of  Engineer- 
ing Mechanics  (footnote  4).  The  notation  used  here  is  that  of  Chapters  40  and  43, 
although  the  equations  are  first  made  dimensionless  before  being  written  into  THIN. 

If  time  (t)  is  scaled  such  that 

r = ^Vl 

P 

and  the  displacements,  strains,  rotations  and  stress  resultants  are  made  dimensionless  by 
the  inti  eduction  of 


(v,  w)  = h(v,  w) 

X = hX/R 

(N0,Ng)  = hD(N0.Ng)/R 
CL  = hKQ/R3 


aTm  =llTn,/(R(!+^> 


(e0(O).e0(O))  = h(e0,e0)/R 

(K0,Kg)  = h(K6)k0)/R2 
(M0,Mg)  = hK(M0,Mg)/R2 

(Pr,  P0)  = Eh2(q,  p)/R2 

a0  = h0/(R2(l+^)) 


where 


D = Eh/(1  -v~) 


K = Eh3/(  12(1-i>2)) 

(Kr,  K0)  = wR:Ar(kr,  k0)/Eh 


h/2 

Tm  = s/  <t-trepw 

•'i  //-» 


e = r|  /*  (T-TREP)fdt 

h 4/2 


it  follows  that  the  governing  equations  summarized  in  footnote  4 become 


^ = w + 30 


€o  = W + V COt  <t> 


V _ 3w  _ 

X “ 'd<t>  ~ v 


- _3X 

K<p  d<j> 


Kn  = X COt  0 


(N0,Ng)  = (e0,§g)  + I2(i0,i0)-  Tr 
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Appendix  C 

THE  TEMPERATURE  INTEGRALS 

In  this  section,  formulae  will  be  derived  for  the  numerical  evaluation  of  the 
temperature  integrals  Tm,  0 which  appear  in  the  constitutive  equations.  It  is  assumed 
that  temperature  data  are  available  at  each  meridional  station  at  five  equally  spaced 
points  through  the  thickness  (see  Figure  4)  not  including  the  end  points  of  the 
interval.  This  latter  restriction  arises  from  the  fact  that  temperatures  are  generally 
computed  using  an  algorithm  based  on  energy  balances  applied  to  regions  of  space.  If 
these  regions  are  chosen  of  equal  thickness,  the  temperatures  are  naturally  obtained  at 
the  points  noted  in  Figure  4.  Alternatively,  should  the  temperatures  be  defined  at 
equally  spaced  points  through  the  thickness  including  the  end  points  of  the  interval, 
the  equations  given  below  can  be  replaced  by  the  well-known  Newton-Cotes  formulae. 

Assuming  now  that  the  temperatures  are  defined  as  shown  in  Figure  4,  we 
begin  by  attempting  to  represent  the  temperature  distribution  by  the  polynomial 

T(j-)  = a0  +alr  + a2F+a3F  +aj4  (f  = f/h, < f < {) 

and  evaluate  the  coefficients  aQ,  . . . a4  so  that 

T,  = T(-0,4)  = aQ  - 0.4a|  +0.42a2  - 0.43a3  + 0.44a4 
T2  =T(-0.2)  = a0  - 0.2a,  +0.22a2  - 0.23a3  +0.24a4 
T3  = T(0)  = aQ 

T4  =T(0.2)  = a0  + 0.2a,  +0.22a2  +0.23a3  +0.24a4 


T5  - T(0.4)  = a0  + 0.4a,  + 0.42a2  + CL43a3  + 0.44a4 


The  coefficients  are  readily  obtained  from  equations  obtained  by  adding  and  subtrac- 
ting the  first  and  last  pairs  of  the  equations  given  above  so  as  to  separate  the 
unknowns  a-,,  a4  and  a,,  a3. 

It  then  follows  that 

a,  = (T , - 8T2  +8T4-Ts)/2.4 

a2  = (-T,  + 16T2  - 30T3  + 16T4  - T5)/0.96 

a3  = (-T,  + 2T2  - 2T4  +T5)/0.096 

a4  = (T,  - 4T2  + 6T3  - 4T4  + T5  )/0.0384 
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Further,  with  the  result  that 


we  obtain 


(275T,  + 100T2  +402T3  + 100T4  + 275TS  )/l  152  - T 


0 = 5(1 +»»); 
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Appendix  D 

THE  STRESS  DISTRIBUTION 

In  this  section,  formulae  will  be  derived  for  the  meridional,  circumferential 
components  of  the  in-planc,  normal  stress  that  are  consistent  with  the  definitions  of  the 
stress  resultants  and  the  sign  convention  shown  in  Figure  1 . 

Let  us  begin  by  defining  the  meridional,  circumferential  components  of  strain  at  an 
arbitrary  point  (</>,  f)  in  the  shell  as  follows 

e6  = e0(o)(0)  - fK0(0) 

Further,  let  the  stress  resultants  Nj,  be  defined  as  integrals  of  the  appropriate 
components  of  normal  stress  through  the  thickness  of  the  shell  with  the  convention 
that 


h/2 

-h/2 

h/2 

(M<p,Me)  = - f (o00,o00)?dr 
-h/2 

Finally,  if  the  state  of  stress  in  the  shell  is  nearly  that  corresponding  to  plane  stress 
(o^  = 0),  then  Hooke’s  Law  takes  the  form 

= - wM)/K  + «(T-TREF) 

e0=(a06  - ^0)/E  + a(T-TREF) 

Alternatively, 

+ «»)/<  l-»2)  - «RT-Tref)/(I-*) 

a00  = E(e0  + ve(())l(\-v2)-  aE(T-TREF)/(  1-r) 

These  equations  can  be  rearranged  in  a number  of  ways  to  suit  one’s  convenience. 

If  the  equations  for  o00,  Oqq  are  substituted  into  the  definitions  of  the  stress 
resultants,  and  we  further  substitute  the  expressions  for  the  strain  components  at  a 
general  point  in  terms  of  the  middle-surface  strains  (e0(O\  and  curvatures  (k0, 

Kq)  we  obtain  the  following  constitutive  equations 
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N0  = D[e/o>  + «/o>-(l+^T11,] 


N8  = D[£(l<0>+1.e/0>-(l«)aTm] 


M0  = K[k^  +vk6  + (l+r)a©3 
M0  = K[xg  + vk^  + (l+r)a©] 


which  are  consistent  with  those  equations  noted  in  Appendix  A.  These  equations  can 
be  rearranged  to  yield  the  following  strain-force  stress  resultant,  curvature-moment 
stress  resultant  equations: 


^(0)  = (N^Ng)/Eh  + aTr 


e0(O)  = Eh  + aT, 


= 12(M^-i^Mg)/Eh3  - a© 
Kg  = 12(Mfl-i,M0)/Eh3  -a© 


Finally,  if  these  expressions  for  the  strain,  curvature  components  are  substituted 
into  Hooke’s  Law,  we  obtain  the  following  expressions  for  the  in-plane  stresses 


°H/E  = VEh  - 1 2rNyEh3  - afT  - Tm  - TR  EF  - f 0)/(  I -k) 


ale /E  = N„  /Eh  - 1 2fM0  /Eh3  - cfT  - Tm  - TR  E F - ?©)/ 1 -►) 


When  expressed  in  dimensionless  variables,  these  equations  become 


<WE  = + Tm  - rS(MF  - e>]/<i-^2)- «(t-tref)/<i-,) 


o#9/E  = R[N#  +Tm  - fh(M„  - 0)]/<l-i’2)-  0<T-TREF)/(l-h) 


f 

§T 


Appendix  E 

ASYMPTOTIC  SOLUTION  FOR  A THIN  SHELL 


In  this  section,  we  consider  the  problem  of  determining  the  deflections,  stresses 
in  ii  clamped  hemispherical  shell  that  is  loaded  by  a uniform  distribution  of  normal 
pressure.  In  particular,  we  will  seek  a solution  as  a sequence  of  solutions,  each  term  in 
the  sequence  multiplied  by  an  ever  increasing  power  of  the  thickness  parameter  (h),  so 
that  (one  expects)  the  accuracy  of  the  solution  at  any  stage  should  increase  as  the 
number  of  terms  in  the  sequence  is  increased.  Although  this  method  of  constructing 
solutions  is  described  in  a number  of  texts,  we  will  refer  here  to  the  formalism 
described  by  Cole.8 


The  problem  we  wish  to  consider  is  described  by  the  equations  presented  in 
Appendix  A where  we  take  pr  = constant  and  delete  all  temperature  terms.  We  wish 
to  obtain  a solution  subject  to  the  boundary  conditions  that 


v,  w,  X = 0 0 = (j)Q 

As  a starting  point,  we  observe  that  the  following  membrane  state 


N0,N„=prR/2 
v = C((1 1 sin  0 


X-M0,m0)q0  = o 


w = ^ - C(0  ] cos  0 


is  an  exact  solution  of  the  governing  equations,  although  it  does  not  satisfy  the 
boundary  conditions.  Physically,  we  expect  that  the  solution  we  seek  should  behave 
like  this  membrane  state  in  regions  of  the  shell  that  are  not  near  the  edge.  What  is 
required,  therefore,  is  a solution  that  is  valid  near  the  edge,  satisfies  the  boundary 
conditions  and  converges  to  the  membrane  state  in  some  sense  for  0 such  that  0<0Q. 

In  this  edge  region,  we  expect  a rapid  variation  with  distance  for  the  dependent 
variables.  In  particular,  we  assume  that  the  dependent  variables  may  vary  significantly 
over  a distance  comparable  with  the  characteristic  length  (v/Rh),  and  define  a 
dimensionless  length  scale  (£),  measured  from  the  edge,  given  by 


0 = 0O  - e£  e2  = h 

Let  us  further  adopt  the  following  set  of  dimensionless  variables  ' 
v = eprR2v/Eh  w = pfR2w/Eh 

(N0,Nfl,Q0)  = prR(N0,N0,eQ) 


u 

Cole,  J D.  Perturbation  Methods  in  Applied  Mathematics , Blaisdell  Publishing  Co.,  1968. 
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<M0’  M0}  = Prh2(^0,  Ue)K\2e2(\-v2)) 

In  these  variables,  we  can  develop  the  governing  equations  in  the  following  form 
(1-p  )(N0,  N^)  = ( l+r>)w  - ( 1,  v )^|-  + l)ev  cot  0 + 0(e3) 

(M0>  M6i)  = d>  ~ ^ cot  <t>  + (1.  «^)e2^|  + 0(e3) 

d^0 

In  V0  = - ~dT  + 6W0-  ) cot  0 + 0(e3 ) 

dN0  _ 

~djr~  e(N0-N0)cot0  + e2Q  + O(e3)  = O 

/s 

uQ  /N  ✓S 

dT"  cQcot0-  (N0+N0)  + O(e3)  = -1 

where  m4  = 12(lV),  and 

cot  0 = cot  0Q  + etj  csc20o  + e2£2  cot  0Q  csc20o  + Ofe3) 
as  appropriate. 

From  the  form  of  these  equations,  it  is  apparent  that  expansions  of  the  type 
v(e,  *)  = v(0)(£)  + ev0 >(£)  + e2  v<2  >(£)  + 0(e3) 
w(e,  £)  = w(0)(£)  + ew(1)(J)  + e2w(2)(£)  + 0(e3) 

N0(e,  £)  = 1/2  + en0(1  >(£)  + e2n0(2)(£)  + 0(e3) 

N0(e’  £>  = n0(O)(£)  + en0(1  >(£)  + e2n0(2)(£)  + 0(e3) 

Q(e,  £)  = q(0)(£)  + eq(1  }(£)  + e2q(2)(£)  + 0(e3) 

M0(e,  = m0(O>(^)  + ™0(1 >(£)  + e2m0(2)(£)  + 0(e3) 

M0(e’  £)  = mC|(O)($)  + em5(1)(J)  + 62m0(2)(f)  + O(e3) 

should  exist.  The  constant  (1/2)  as  the  leading  term  in  the  expansion  for  N 
anticipates  the  requirement  for  zero  order  matching.  ^ 
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The  equations  governing  the  zeroth,  first  and  higher  order  terms  in  the  above 
expansions  are  obtained  by  substituting  these  expansions  into  the  governing  equations 
and  taking  successive  limits  of  the  governing  equations  as  e ->  0 holding  £-fixed.  In 
particular,  we  find  that 


dvf0) 

djt  = n+n)w(0)  - (lV)/2 

=(!+„)»«>>  - ^H<0> 

(0)  _ 

df  ne  ~ 

nyo)  = d2w(0) 

o 

w-o- 

£ 

u 

o 

E 

df2 

govern  the  zeroth  order  solution,  and 

(l-i'2l(n0,1).nj,"l)>(l+^®">-  ’+(■.,  l)v<0»  cot 


«"V  ’>=0,^1 


(i) 


dw 


(0) 


mV' 


dn0(1) 


dm. 


(i) 


- 1)-^  cot  0O 

df  +(m0(O)  - m{?(O))cot0o 
df  ~ (U~~  n^(0))  cot  0Q 


l ) 

df  q(0)  cot  0O  - (nf l)  + n0(1))  = 0 


govern  the  first  order  correction,  and 


(1  Xn0(2),no(2>)  = n ^)w( 2 > - ( 1 , *>)£  ‘ + ( v;  1 )^° J > esc2 0O  + (n,  1 )v< 1 > cot  0 


'df 


fm(J(2),  mf2)) 


,,  ./d2w(2)  dv<0)\  / dw(0)  , dw(!)  \ 

’l'  (d^2  + df  /"^’^df  csc^o  + d|-  cot  *o) 
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v2) 


dm  j ^ ^ 

- "dT  +Hm0(O)-me(O))csc20o  +(m0(,)-m0(]))cot0o 
q(0)  +(l/2-n0(O))  £ csc20o  + (n^1  )-n0(l })  cot  0Q 


dq(2) 

d£  " ($q(0)  esc 20o  +q(1)  cot0fl)-  (n0(2)  + n0(2))  = O 

govern  the  second  order  correction.  In  the  paragraphs  to  follow,  the  solution  to  the 
a ove  system  ot  equations  will  be  described  which  matches  the  membrane  state  in  the 
region  of  the  shell  that  is  not  near  the  edge.  It  will  be  assumed  that  the  reader  is 
familiar  with  the  matching  process  so  that  little  of  the  motivation  needs  to  be 
described  and  only  the  results  of  the  process  will  be  presented.  Furthermore,  only  the 
case  0O  - 90  deg  will  be  described. 


Thus,  in  anticipation  of  the  requirement  that  all  solutions  be  bounded 
exponentially  as  £ ->  °°  we  take  the  zeroth  order  solution  in  the  form 


w(0)  + e-x(A0  sin  x + BQ  cos  x) 


x = m£  A/ 2 


q(0)  --e  X[(A0  +B0)  cosx  + (A0  - B0)  sin  xj/nv^ 
v(0)=C0(i)-f(l-f,)q(0)  n0(O)  = u/2  + w(0) 


The  solution  of  the  equations  governing  the  first  order  correction  can  be  shown 
to  vanish  identically  for  the  case  0Q  = 90  deg.  This  follows  from  the  anticipated 
observation  that  as  the  governing  equations  become  homogeneous  for  0n  = 90  deg,  the 
on  y solution  which  satisfies  homogeneous  boundary  conditions  is  the  trivial  solution. 

Thus,  as  we  seek  more  than  the  first  term  in  the  expansion,  we  must  turn  our 
attention  to  the  second  order  equations.  In  doing  so,  we  find  that 


w ] (l+i>)C2(l)  C0(,)£- 3£q(0)/4  + e"x[£2(AQ  sin  x + BQ  cos  x)/4  + 

+ A2  sin  x + B2  cos  x] 
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?(2)  = D,(i)  - 2(l+I>)C2(i)£  - C0(i)£2/2  + ~ e"x  sin  x(— 


m 


A2  , Ao  + Bo 

+ 7 5 h 


4m 


j 


V B„-A0  a / 

- i + ~ — 7-  X-  +COSX  - 

2 nr  2m3  / \ 


n/2,  = ™/2>  + w<2>  + £v<0) 


A2+B 2 , B0  - A0  Aox  A0  + B0  , 


m 4m3  2m3 


^x2) 
2m3  / 


The  constants  Afl,  B0.  CQ(i\  A,,  B-,,  C2<*\  D2(i)  of  the  edge  (or  inner)  solution,  and 

( of  the  membrane  (or  outer)  solution  will  now  be  chosen  so  as  to  satisfy  either 
the  boundary  conditions  or  the  matching  requirements. 

The  boundary  conditions  at  0 = 90  deg  (£  = 0)  require  that 
w(n)(0)  = w(1  *(0)  = w(2)(0)  = . . . = 0 


v(0)(O>  = v(1  )(0)  = v(2)(0)  = . . . = 0 


dw(0) 

<1* 


dw(2 } 

d* 


(0)  + v 


(0) 


(0)  = . . . = 0 


On  applying  these  boundary  conditions  to  the  zeroth  order  solution,  we  find 


A0  = B0  =-(l-i»)/2  c0(i)  =-(\-v2)ls/2m 

Now.  before  applying  the  boundary  conditions  to  the  second  order  corrections,  it  is 
convenient  to  obtain  the  requirements  of  the  matching  process. 

Following  Cole  (see  footnote  8),  we  first  define  an  intermediate  length  scale  £ 
where  4 


_0O-0 

7?(e)  - &lv(e) 


and  77(e)  is  such  that 
e-*0 
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For  example,  we  could  chooser?  = e’^2.  In  what  follows,  we  will  require  that  the 
edge  (or  inner)  solution,  expressed  in  the  intermediate  variable,  converge  to  the 
membrane  solution,  also  expressed  in  the  intermediate  variable,  for  all  values  of  £ 
(fixed)  as  e-K). 


•v 


On  applying  this  matching  condition  to  N^,  we  find  that  n^'^-O  thus 
verifying  the  fact  that  the  first  order  equations  are  homogeneous,  and 


L 


y2>=o 


e-*0,  - FIXED 


so  that 

C2<#-0 

On  applying  the  matching  condition  to  v,  we  find 
C(o)  =eprR2C0(i)/Eh 

for  zeroth  order  matching.  The  remaining  constants  are  found  on  applying  the 
boundary  conditions  to  the  second  order  correction. 


dw1'  • 

The  boundary  conditions  on  w(2\  are  satisfied  by  taking 


(2) 


A9  = -(l-y)(l+4y)/4m- 


B2  =0 


The  boundary  condition  on  v*2*  is  satisfied  by  taking 


D2(i)  = -m(l+4r0/48v/2 

This  completes  the  solution  as  it  can  be  shown  that  the  remaining  matching  conditions 
are  also  satisfied. 

As  a final  remark,  it  is  convenient  to  obtain  the  above  solution  in  a fonn 
which  is  uniformly  valid  over  the  entire  range  of  opening  angle  rather  than  have 
solutions  whose  range  of  validity  was  restricted.  Again  Following  Cole  (footnote  8), 
such  uniformly  valid  solutions  can  be  constructed  by  subtracting  the  common  part 
from  the  sum  of  the  inner  and  outer  solutions.  The  common  part  of  the  solution  are 
those  terms  in  both  the  inner  and  outer  solutions  which  match  identically  in  the 
matching  process.  For  the  example  studied  here  (0O  = 90  deg),  it  can  be  shown  that 
uniformly  valid,  two-term  expansions  can  be  constructed  in  the  following  form 
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Appendix  F 
PROGRAM  LISTING 


EFFECTS 


QB  EQUALS  Q 
E PPn  EQUALS 
EQUALS 
EQUALS 
EQUALS 
E QUALS 
EQUALS 
EQUALS 
EQUALS 
EQUALS  W 
EQUALS  V 


EPTR 

FKPP 

FKTn 

FNPP 

FNTP 

fmpb 

FMTD 
WA 
VA 


THIN, SPHERICAL  SHELL  WITH  SURFACE  LOAD , TEMPERA TURE 

boundary  T°rn  1°"  5L*MPEC  BOUNOm  *T  PHI=PO.  TO  1 FOR  SIMPLY  SUPPORTED 
BOUNDARY,  TO  2 FCR  ROLLER-SKATE  BOUNDARY. 

glgj  rTorcrr°  °Tn°-  ° ISPLAC EMENTS  ONLY. TO  1 FOR  BOTH  DISPLACEMENTS 
AND  STRESSES,  TO  2 FOR  STRESSES  ONLY. 

HB  EQUALS  H/R. 

CKR  EQUALS  KR 

CKP  EQUALS  KP 

TMB  EQUALS  TEMPERATURE  MEAN  BAR 

THB  EQUALS  TEMPFRATURE  FTPST  MOMENT  BAR 

PB  EQUALS  P BAR  (TANGENTIAL  PRESSURE) 

QPB  EQUALS  Q BAR  (TRANSVERSE  PRESSURE! 

WP  EQUALS  W BAR 

WBD  EQUALS  U BAR  DOT  (VELOCITY) 

VB  EQUALS  V BAR 

VBD  EQUALS  V BAR  DOT  (VELOCITY) 

CHIB  EQUALS  CHI  BAR 
BAR 

FPSILON  PHT  BAR 
EPSILON  THETA  BAR 
KAPPA  PHT  BAR 
KAPPA  THETA  BAR 
N PHI  BAR 
N THETA  BAR 
M PHI  BAR 
M THETA  BAR 
BAR  NEXT  TO  APEX 
BAR  AT  A STATTON  I=NV. 

CHOOSE  N CRCATER  THAN  OR  EQUAL  TO  ( i» 2 .PO/SQRT (HB M /2 
CHOOSE  NV  ABOUT  N/2. 

DIMENSION  PB  ( 5 C ) « 8PB I 50)  »TNB( 50) , TMB I 501 
TMP(  50,  5)»ST(5BI5I 
CT1(50),CT2(50) 

HB(  50 ) . WBD( 50  I, VBI 50) »VB 01 50 1 .CHIB  I 50  I 
03(50), FNTBI 50 ) ,FNPBI  50 1 , FMPB I 501 »FNTBI50) 

WA(  3500  ) . VA(  3500 ) 

geometric, material  and  control  data  entered  mere. 

RE  AD (5,24)  PO , HB , PP , CK R »CKP 
READ ( 5 « 24  | TPET.ALFA 
READ (5,25)  NBC, N, NV. MMA X, IDEC,  IPRNT 
FORMAT ( ) 

FORMAT (6110) 

L = N-1 
AN  = N 

DP  = 2,OPC/( 2.0*AN-1.0I 
0T=0.5*DP*DP 
S1R=1.0-0.5*CKR 
SlP=1.0-0. 5*CKP 
S2R=1 .C.G  .5  *CKR 
S2P=1.0*0.5*CKP 

PRESSURE , TEMPERATURE  DATA  ENTERED  HERE. 

TEMPERATURE  INTEGRALS  EVALUATED  ASSUMING  TE*»ER  ATURES  DEFItCD  AT 

FIVE  EQUALLY  SPACFD  POINTS  THRU  THICKNESS,  NOT  INCLUDIN6  SURFACE  POINTS. 

J=l,5  REFER  TO  POINTS  NEAREST  INNER.OUTER  SURFACE  RESPECTIVELY. 

READ ( 5 , 24  ) ( QP8 ( I J » PB ( T ) .1  = 1 ,N ) 

DO  20  1=1, N 

RE  AD  ( 5 , 24  ) ( TUP  ( I , J ) , Jr  1 ,5  ) 

, J * = * l*  °*PPI  .ALFA  •(  - TREFH  2 75.0  *TMP  1 1. 1 )» 100.0*  TMPI 1 ,2 1 .402.0 

1«TMP(I,j)»ioC.O.TMP(I,4)»275.0*  TMP I 1,5)  1/1152 .01 /MB 


DIMENSION 
DIMENSION 
DIME  NSION 
DIMENSION 
DIME  NSION 


24 

25 


45 


ammmm 


■■mmm 


m*rnmm 
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THBm  = II.  O.PR ) *ALFA • ! - U.O*TMPII*  11-2.0 »TMPII*2|  »2.  0*TMPI  1.41 
1»11  .C»TMP II , 5 I )/ 19.G*HC*HB1 

20  C ONT INUr 

WRITE  16.21)  PC.PR.CKR.CKP 
WRITE  I E « 22  1 NBC.N.NV.MMAX.HB 

21  F 0 R M A T I 1 H •,P0=,»E11.5«2X.,PR=*«Ell.5»2X.,KR=,»E11.5.2X. 
1*KP=’.F11. 51 

22  F0RKATI1H  . *NB C= * . I I . 2X . »N= • . 12 . 2X.  'N  V = ' . 12 • 2 X . ' MMAX= ' 1 1 4 • 2 X • 

1 •HB=*.ril.5) 

WRITCI6.23)  TREE. ALFA 

23  F0RMATI1H  . * TREF= • . Ell . 5. 2X • • AL FA=  • . El 1. 5 1 
WRITE  I S • 4 C I 

40  FORM  ATI //, 5X. 'TEMPERATURE  DISTRIBUTION • /) 

WRITHE.  35  1 / 

DO  42  1 = 1. N / 

WRITE ( G . 4 1 1 I.  ITHPIT. J)i  J=1.5I 

42  C ONT INUr 

41  Format i ih  , 12, 2x .5 < eii.5 . 2x i i 
WRITE  16.431 

43  F0RMAT(//,5X . 'SURFACE  LOADING  ANO  TEMPERATURE  INTEGRALS • / I 
WRITE  10.44  1 

4 4 FORMAT  I2X .*1 ». 9X . • GPB • . 9 X . • PB  • . 9X . • TMB  • • 9X » ‘THB'/Zl 

DO  4 G 1 = 1. N 

WRITE  IE. 4 0 I I. GPB 1 1 1 . PBl T 1 . TMBl 1 1 . THB 1 1 1 

45  F0RMATI1H  .1 2 . 2 X. 41 E 11 . 5. 2 X I) 

46  CONTINUE 

SET  COTANGENT  FUNCTIONS 

CT2  II l=r0SIDP/2.CI/STNlOP/2.O> 

DO  1 1 = 2. N 
FI=  I 

CTllII=roS(IFI-l.  ai*0P1/STNI  IFI-l.  01  *DP  I 
CT2  (II=COSI I PI -0.5  I »0P l/SINI I FI-0.5  I *0P» 

1 CONTINUE 

SET  TIME  TO  ZERO 
M=  1 

SET  INITIAL  DISPLACEMENT. VELOCITIES  TO  ZERO 
DO  2 1=1, N 
WBII 1=0.0 
WBD 111=0.0 
VP  1 1) =0.0 
VDD 111=0.0 

2 CONTINUE 
W A 1 1 I = C • 0 
VA 111=0.0 

SET  APEX  SYMMFTRT.  TANCENTIAL  DISPLACEMENT  BOUNDARY  CONDITIONS. 
CHI8I  11=0.0 
GBI1  1 = 0.0 
VP  I N * 1 1=0. a 

DISPLACEMENTS  COMPLETE  FOR  FIRST  TIME  STEP 

COMPUTE  STRAINS. ROTATIONS. STRESS  RESULTANTS  FOR  TIME  STEP 

3 DO  4 1 = 2. N 

CHIBII1=IUBII l-WBII-1 I l/DP-VBl 1 1 

4 CONTINUE 
DC  5 1=1. L 

F PPB  = UP 1 1 1 VBI  I* 1 l-VBI II J ZOP 
EPTBrWHI  1*0. 5*  I VBIT*1I*  VBI  II )*CT2I  II 
FKP3  = I EHIBII  ♦U-CHI9I  T1  1/DP 
FKTB=C.5*ICHID  11*1 1 *CHTB  III  1 *CT2m 
FNP3 II )=EPP3.PR*EPTB-TMBI I I 
FNTBIIl=EPTB»PR»EPPB-TMBIII 
TMP3 II 1=FKPB*PR*EKTB*  THB I I 1 
FMTBIII=FKTE*PR»rKPB*THBITI 

5 CONTINUE 


46 


..  aii 


t 


ICO 


?nc 


3D0 


)-  THBI  N|*  DPI  / 


1 I 


NWC  TP  5785 

SETT  CONDITIONS  AT  I=N. 

FPP3  = W'M  N)  «•(  VB(  N.  1 )-V9(  Nl  I /DP 
F P TR : W° ( N I *0 . S • ( V3 ( N ♦ 1 )♦  V3| N I I • C T 2 1 N J 

FNP3lMI=FPPB.P9.EPTB-TKBINI 

FNTR  ( N I = F P T G ♦ P R » ("  p p p - y Rg  j pg  j 
iriNBC-ll  IOC, 200,200 
F KPR r-D . c *c HI 3 ( N ) /□ p 
FPP9(N)=FKPB*THB(  Nl 
FPTCINl=PR.rKP8»THB(N) 

GO  TO  300 
FKPB f N I =0.0 

C HIP  { N ♦ 1 , r | CHI  9 ( ,N  ) •(  1.0-0. 5.P9. DP. CT2»  Mil- 
l(l.C*D.C*PR«DP»CT2lN)  | DP  * CT  2 f N I 

FKPP,  = (CHI3(N*1)-CHIB(  N I )/op 
FKTRrD.5»(CHIB ( N ♦ 1 I ♦ C H T B ( N I l*CT2(NI 
FPT3(N)=FKTB«-PR«FKPB*TM3’N| 

00  E I=?,n 

i-rMTR  a-nflcnrir8* 1-1  ,,/dp*d*5**fnpb iii-fhtbi n .fiipsi  i- 

E C ONTINUr- 

C0MS?T7  STEPTED  F°R  ™IS  TIHE  STEP 

J*G*5*  Jr-VJS  I*  -- N f B f T TI FN  P ^ |Sr  - 1 1*- ? ®T*  J ^ ‘ 1 f T T f * “FNP0  ' ^ 1 » >'°P 

^CCNTINUF3'1  ,/12*0  ,/(  1,0“P***P*M 
oo  3 1 = 1, l 

2 l/l l.r-rR.PR| I tl 1/ DP  10*5*1 QB 11*1 1 '06  1 1 ) > • c T2 1 1|  1/12. o 

8 CONTINUr 

IFIMBC-ll  32c.320.310 
310  Qn.  (NH  I =-QBC  Nl 

?l/(I.O-PR.PR)  , ’ DP  *°*5*»8B»N*1  >*QB(NI  I.CT2INI  1/12.  D 

320  P = M»1 

COMPUTC  RISPLACFHFN  TS  TOR 
WRI1IZWPI1 I *0  T»  WB  Ot 1 ) 

WA  ( !', I =WE  ( 1 1 
00  9 1 = 2, L 

wp 1 II -W° (I| ♦DT«W8D ( I I 
VP (I I = VE (I l»DT*VBD( Tl 
9 CONTINUr 

VBIN)  = VE(NJOT«VBD( Nl 

VB  ( N ♦ 1 ) r -V9  ( N I 

VA  (Ml  = V?  ( N VI 
IFINPC-IJ  3M0, !i(0,33C 
WP(N)  = WP(N|.OT»WDO(  Nl 
IF(P-PMAX)  3.3,10 
IFIIPRNT-1 | U,ntli, 

WRI  IF (6 . 1G I 

FORMAT.//, 2X.. DISPLACEMENTS  VS, 

Oc  13  J = 1 0, MM A X » IDEC 

19  ^niTE:, 6,121  J,HA*  J»iJ«VAIJ  1 

I,  cKl?"  

IF(IPRNT-1|  15.1K,  m 
19  WRITE (6,30) 

!n  "!0lUTI°"  ” »I»EKIO«XK  ..RI.8LESVI 


NEXT  WHOLE  TIPE  STEP 


330 

3NC 

10 

11 

IE 


ITERATION  INDEX*/! 


47 


NWC  TP  5785 


DO  12  1=1, N 

32  CONTINUE 
WRITE  IE . 3 4} 

" Ei"s:irj;r-Mo,,“— • mew**.,, 

DO  37  1 = 1,  N 
DO  36  J=l,5 
ft  J = J 

ZETB=-D.E*0.2*AJ 

36  CONTINUE 

WRITEJE.mi  I,(SIG8IJI,J=l.5| 

37  CONTINUE 
15  STOP 

END 
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